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Abstract. We review our recently developed electronic structure calculation 
methods for the dynamics of large-scale solids, liquids or soft materials with an 
efficient algorithm of linear algebra. The Hamiltonian of electronic structures is 
based on the 'atomic superposition and electron derealization molecular orbitals 
theory' (ASED), using the Mulliken charge density. Very crucial algorithms of 
the linear algebra are the generalized Lanczos method and the generalized shifted 
COCG (conjugate orthogonal conjugate gradient) method for general eigen-value 
problems. The shift equation and the seed switching method are essential for 
the shifted COCG method and reduce the computational cost much. We, then, 
present some applications to electronic structure calculations with MD simulation, 
such as the fracture propagation in nano-scale Si crystals and the proton transfer 
in water. 



PACS numbers: 61.46.-w,71.15.Dx,71.27.+a 
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1. Introduction 

The electronic structure calculations in nano-scale structures have attracted 
considerable attention. It is much desirable to achieve molecular dynamics (MD) 
simulations with quantum mechanical electronic structure calculation for systems of 
several hundred thousands atoms with a few hundred pico-seconds (or longer time) 
process, since the atomic structure is very essential to the electronic structure and 
vice versa. Requirement for higher accuracy of the total electronic structure energy 
in nano-scale systems and that for larger system size contradict each other. One of 
possible choices for satisfying these requirements could be the LDA (local density 
approximation) calculation with massive-parallel computation. Even in that case, a 
long-time process for MD simulations is another difficulty. 

In the present paper, we propose one possible way to pursue a quantum 
mechanical MD simulation of large system size and long time process at the expense 
of lower accuracy of the electronic structure energy. The first-principles tight-binding 
method is one of promising techniques and we review the "atom-superposition and 
electron- derealization molecular- orbitals theory" with self-consistent charge density in 
Section \2. 21 Then, based on the tight-binding formulation, we develop novel algorithm 
of linear algebra for large scale systems with the non-orthogonal local orbitals in 
Section [3] Section 2] is devoted to present applications to nano-scale systems, the 
fracture propagation and surface reconstruction in Si single crystals and the proton 
transfer in water. Section [5] is conclusions. 

2. First-principles tight-binding methods and charge self-consistency 

2.1. Several first-principles tight-binding methods 

We have several semi-empirical tight-binding formalisms [U [2] whose results are 
comparable to those of the first-principles LDA calculation. We have also several 
theories of first principles tight-binding Hamiltonian, e.g. the tight-binding LMTO 
method, [3J the LDA tight-binding method, [3] the first principles tight-binding 
method of the Naval Research Laboratory (NRL) . [5] We explain in the next subsection 
the tight-binding Hamiltonian of the 'atom-superposition and electron-dclocalization 
molecular-orbital' theory. [6] 

2.2. Atom-superposition and electron- delocalization molecular- orbital theory 

The atom- superposition and electron- delocalization molecular- orbital (ASED) theory is 
based on a charge-density partitioning method, where the charge density is determined 
by the Mulliken method of partitioning the overlap charge density (the Mullikcn 
charge) of assumed atomic orbitals. 

Diagonal elements of the Hamiltonian matrix in the extended Hiickel theory would 
be given as the valence orbital ionization potential with a quadratic formula of the 
Mulliken charge. In the present form of the ASED theory, a diagonal element of 
the Hamiltonian Hi a .i a , where i and a denote an atom and its orbital, is set to 
equal to an atomic energy level (ionization energy) ei a . An off-diagonal element of 
the Hamiltonian is determined by the bond formation and an extended Huckel-likc 
delocalization energy may be generally a good approximation; [7] 

tt zs(„ a. D -\ {Hia,ia + Hjfjjp} 
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where a coefficient K(a, /3; -R«) stands for distortion of wavefunctions depending upon 
orbitals and a distance [5] and S ia _jfj is the overlap integral. 

Another contribution to the total energy is the ion-core repulsive energy, which 
is defined, in the tight-binding density functional (DFT) formalism [5], to be the 
difference between LDA energy and the tight-binding band structure energy. 

The charge transfer effects is crucial within the charge self-consistent (CSC) 
theory and formulated by the second order perturbation theory, 10; which gives a 
form of the extended Hubbard-type Coulomb interaction as 

H? a % = -SfaJ/J J2^ k + ^)Aq k (2) 
k 

where Aqk is the deviation of the (Mulliken-type) charge of an atom k and ^ik is the 
(diagonal or off-diagonal) Coulomb interaction between atoms i and k or the chemical 
hardness. 

We are applying the CSC tight-binding method to water, molecular liquids and 
soft materials of nano-scale size (~ 3000 atoms). A MD simulation of a proton transfer 
in water will be presented preliminarily in this paper. 

3. Linear algebra and generalized Krylov subspace method for 
generalized eigen-value problem 

3.1. Linear equation and Green's function [llj 

In order to study electronic structures in nano-scale materials, it is much desirable 
to construct the order-A algorithm whose computational load, e.g. cpu time and 
memory size, is proportional to the system size or the number of atoms. 

Problems on electronic properties in materials starts usually from obtaining eigen- 
energies and eigen-functions in materials or solving linear equations to find |a;")) for 
a given state \ j) and a given energy E\ 

(E-H)\x^) = \j}, (3) 

where % is a one-electron (Hermitian) Hamiltonian operator. Once we solve Eq.Q 
and find a vector \x^'), the matrix elements of the Green's function G(z) = (z — 7i) _1 
can be evaluated very easily as 

Gij (E) = (i\(E + i5- ny'lj) = (i\x^), (4) 

where 5 is an infinitesimally small (positive) number. 

3.2. Generalized Lanczos method |12) 

In case of the orthogonal basis set, the Lanczos process and the shifted conjugate 
orthogonal conjugate gradient (COCG) method [13j HH [15] can be constructed in the 
framework of the Krylov subspace generated by % and x^ '. 

When the basis set {(pi} is non-orthogonal, the overlap of wavefunctions \tp a ) 
(IV-'q) = J2i l^i) - "!^) is define by the "S'-product" as 

<vg^) = («(<*U (/9) ) ^uW'ufstj (5) 

ij 

and Sij = {4>i\4>j) is an element of the overlap matrix. 
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With a non-symmetric (non-Hermitian) matrix H — S~ 1 'H, we can get the three- 
term recursive equation or the generalized Lanczos process as 

Hu n = a n u n + b n+1 u n+1 + b n u n ~ x (6) 

for n = 0, 1, • • ■ and bo — with the "S-orthogonality" as 

(u n+1 ,u m ) = for m = n,n-l,n-2, ,0. (7) 

We can set b n 's (n = 1, 2, • ■ •) to be positive numbers. The generalized Krylov subspace 
is spanned by the vectors 7i m x^ — (S^ 1 H) m x^ and defined as 

Kn(H, z (0) ) - Bpan{s<°>, Hx®\ #V°), . . . , (8) 

Within the set of basis {it , it 1 , u 2 , ■ ■ •} or the resultant 'generalized' Krylov 
subspace, the 'Hamiltonian matrix H' is tri-diagonalized on the basis of \tp n ) = 
J2i |^ > i} u i™' ) - Therefore, one can calculate the Green's function in this 'generalized' 
Krylov subspace and also achieve the spectral decomposition (the subspace 
diagonalization) in the 'generalized' Krylov subspace, which we call the generalized 
Lanczos method. 

3.3. Generalized shifted COCG method pl[T2] 

Even for non-orthogonal basis set, we can generalize the COCG procedure, similar to 
the original one, |15] for the Hamiltonian % and the overlap matrix S. The linear 
equation of the 'seed' energy (an arbitrary chosen energy) a s can be written as 

{S- l n + a s l)x {i) = S^b. (9) 

An approximate solution at the n-th iteration, a searching direction for the solution at 
the next iteration step, and a residual vector are written as x n , p n and r n , respectively. 
We can get a recursive equations under an appropriate initial conditions; 

Pn = r' n + P n -XPn-l, (10) 
X n+ i — X n + a n p n , (11) 

r n +i = r n - a n (H + a s S)p n , (12) 

with r' n+1 = S-V+i, On = (p „,(g+ ( r. ) ) p„) and /3„ = '^ff' . 

For non-orthogonal basis set, we can generalize the shift equation. The linear 
equation of a 'shift' energy a is as follows; 

(S^H + al)x {t) = 5 _1 6. (13) 

Then the entire information is delivered to the coefficients a° and j3° at every 
shift energy a by a scalar constant tt^, lf which is generated by the recursive equation 

<+l = {1 + C*„(<7 - «7.)}< ~ — On« ~ <_i). (W) 

OCn-1 

The most important issue of the shift equation is the theorem of collinear residual: |14j 
< = ^r n . (15) 



Then, once the solution is obtained at the seed energy by Eqs. (IT0|) ^ (|T2|) . solutions at 
arbitrary shift energies can be evaluated by using only scalar-sealer and scalar- vector 
multiplications. 
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The only remaining problem is the calculation of S r n . Fortunately, the overlap 
matrix S is positive-definite and nearly equal to the unit matrix. Therefore, we can 
solve this equation efficiently by the CG method, e.g. solve inversely r n = Sr' n for a 
given tv 

The choice of the seed energy a s is not unique. If one would choose a seed energy 
in an energy range of faster convergence, spectra at majority energy points have not 
be converged yet after achieving the convergence at the seed energy. In that case, one 
might restart the calculation with new seed energy from the beginning. 

The seed-switching is very efficient technique to avoid waste of the already finished 
part of constructing the Krylov subspace. [TTl [18] One can choose a new seed energy 
cr" ew and can continue the calculation without discarding the information of the 
previous calculation with the old seed energy a s . 




1e-04 1 • 1 1 ' 1 
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Figure 1. The local density of states of fee gold of Au 256 atoms by the tb-LDA 
with non-orthogonal basis set. There are three lines in the figure and they are 
almost overlapped. The chemical potential locates at +0.2947 Ry. 



3.4- Example 

We have explained the generalized Lanczos method and the generalized shifted COCG 
method. Now we present the calculated results of the local density of states (LDOS) 
by using the tight-binding Hamiltonian with non-orthogonalized basis set given by 
the NRL Hamiltonian. J5j The system is of gold 256 atoms in a supercell of an fee 
structure with s-, p-, and d-orbitals under a periodic boundary condition. Figured] 
shows the LDOS of d-orbitals, where the vertical axis is in logarithmic scale to see 
the behavior of tails. The three calculated results, the generalized Lanczos method, 
the generalized shifted COCG method and the exact calculation, are almost identical 
with each other. 

The cpu times by a standard work-station of a single processor are 3.92 s, 21.20 s 
and 57.6 s for the generalized Lanczos method, the generalized shifted COCG and 
the exact calculation, respectively. The algorithm of the present developed algebraic 
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methods is suitable for parallel computation, on the other hands, the cpu time grows 
as a cubic power of the matrix size in the exact calculation. Therefore, the present 
methods are promising in large scale electronic structure calculations. 

4. Application to nano-scale systems 

4-1- Scale-size dependence of recently developed linear algebraic methods 

Much attention has been paid to nano-scale systems and the first principles 
molecular dynamics simulation has become more important in material sciences and 
engineering. We have developed a set of computational methods of electronic structure 
calculations, i.e. the generalized Wannier state methods, [2U1 US 1221 123] the Krylov 
subspace method, [24] and the shifted COCG method for nano-scale systems. [15] 
In the preceding sections, we have explained recent development on the ASED 
and generalized shifted COCG method. The tight-binding MD simulation method 
with being-developed linear algebraic algorithm is useful to achieve calculations in 
extra-large systems whose computational cost increases linearly proportional to the 
system size as shown in Fig|2j where tested examples are metallic, insulating and 
semiconducting dense materials. 



Figure 2. The computational time as a function of the number of atoms (N). 
12 1 1 1221 1231 The time was measured in a older standard workstation for metallic 
(fee Cu and liquid C) and insulating (bulk Si) systems with up to 11,315,021 
atoms, by the conventional eigenstate calculation (EIG) and by our methods 
for large systems; KR-SD (Krylov subspace-diagonalization), WS-VR (Wannier- 
state variational) and WS-PT (Wannier-state perturbation) methods. The CPU 
time increases linearly with increasing the number of atoms. See the original 
papers 21 22 23] for the details of parallel computation. 



4-2. Fracture propagation and surface reconstruction in silicon crystal 

Here, we will show a calculation results of fracture propagation in a bulk Si single 
crystal with the semi-empirical tight-binding Hamiltonian by Kwon et al 1 , based on 
the orthogonal basis set. The cleavage propagation speed is estimated, in the present 
model, as w pr0 p — 2 nm/ps = 2 km/s, which is comparable to the velocity of the 
Rayleigh wave ur, = 4.5 km/s ~ 4.5 nm/ps. [21] [22] 
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Figure 3. Fracture propagation on a (111) plane of a Si single crystal, (a) 
Appearance of Pandy structure without steps, (b) Pandy structure with a step. 
The numbers show the number of atoms in a ring. The symbols (+) or (— ) refer 
to the electron charge, excess or deficiency. 



The easy-propagating plane of fracture in Si is known widely to be that of (110) 
or (111) and we studied the phenomena of fracture propagation in 14 nm scale Si 
crystals. [22] In case of fracture on the (111) plane, the the Pandy structure [25] 
of (lll)-(2 x 1) surface reconstruction appears with several steps. The direction 
or the structure of fracture propagation planes is not explained by the stabilization 
energy of the final stable surfaces structure because a transient surface structure is 
first formed immediately after the bond-breaking and, after appearance of the ideal 
surface, the surface reconstruction happens to form the final reconstructed surface. In 
the process of fracture propagation and surface reconstruction, the two different kinds 
of energy loss and gain compete with each other, e.g. the energy competition between 
an electronic energy gain of breaking bonds and an elastic energy loss of distorting 
the lattice. This kind competition may be possible only in nano-scale systems and, 
in a MD simulation, we should prepare larger systems, e.g. a few tens nm length. 
Then, even a fracture propagation starts on a (001) plane, the plane of the fracture 
propagation changes to (111) and (110) planes. Figure [3] shows examples of the surface 
reconstruction after formation of surfaces. 



4-3. Proton transfer in water 

The correlation energy of water dimer, [26] hydrogen-bonded network [27] and 
protonated water [28] are all long-investigated problems of quantum chemistry. We 
have been applying the ASED method to water and proton transfer in water to 
investigate the applicability of the method. 

The bond angle H-O-H (an experimental value 104.5°) is determined sensitively 
by the energy levels (binding energy and bond angle dependence) IB2 and 3Ai of 
H2O, which could be done by adjusting the ionic energy levels £2s and e2 P of oxygen 
2s and 2p orbitals and their STO parameters. After that, we could obtain the correct 
order of molecular energy levels of (2Ai), (IB2), (3Ai), (lBi), (4Ai), (2B2) and an 
approximate cohesive energy of a molecule (10.97 eV, Cf. the experimental value of 
10.08 eV). The bond length can be adjusted by slight rescaling the two-body repulsive 
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Figure 4. Snap shot of the proton transfer process. A proton transfers from a 
H2O molecule A to that of B and another proton will transfer from the oxonium 
ion H30 + to that of C. The inset shows the definition of angles a and /3 for a 
stable configuration of a water dimer. 

term and the resultant bond length is 3.168Aand the CSC scheme works efficiently 
(Cf. 2.797Awithout CSC scheme and 2.967Aby Gaussian B3LYP/6-31G(d,p)). Once 
we include H 2p state in the calculation, the molecular dimer configuration a — 7°, 
f3 = 57° can be obtained, which is in good agreement with the results of Gaussian 
calculation (a = 6° and ft = 57°. See the inset in Fig.@}. 

Then we simulate the proton transfer in a water system IOOH2O + HaO + and 
observe the proton transfer process to happen in a recombination of H30 + + H2O => 
H2O + H30 + (hopping of a proton) but not the transfer of an oxonium ion H30 + itself 
and the jumping time is of the order of 1 ps. Furthermore, we observe the following 
elementary process within an order of 0.1 ps; 

(1) A H30 + ion tries to find a proper one in near-neighboring H2O molecules, 

(2) a H30 + ion forms a weak bond with one of H2O (forming H30 + +H20), and 

(3) a proton transfers to that H 2 0, forming H20+H 3 + . 
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(4) If this H 2 is proper one, the new H 2 molecule (the former H 3 + ion) rotates 
(to form a stable water dimer 2H 2 0) and change to a stable configuration and, at the 
same time, 

(5) a new H3O 4 " ejects a proton H + to a next H 2 0. 

From these observation, we may conclude that the ASED framework can work 
properly in a wide variety of materials. The present system size of water is still 
too small to use the generalized Lanczos method and the generalized shifted COCG 
method. Even so, we have already assured these methods for working in the MD 
simulation in the present system. 



5. Conclusions 



We have reviewed our recently developed methods of ASED based on the generalized 
Hiickel approximation and the novel linear algebraic algorithm for large-scale quantum 
molecular dynamics simulation. The crucial point is that the novel linear algebraic 
algorithm is applicable to systems of several thousand atoms or more with the same 
accuracy of the exact calculation. Then we presented examples of the applications to 
the fracture propagation and surface reconstruction in Si and the proton transfer in 
water. The water is very sensitive and difficult systems to simulate and the successful 
results suggest that our proposed simulation method is promising in a wide variety 
of systems, dense or dilute, solids, liquids or soft materials, and metals or insulators. 
The present MD simulation method was applied also to the problems of the formation 
of helical multishell Au nanowire |29j and the diamond-graphite transformation under 
applied stress [5U] . 
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